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Abstract 

Staggered fermions with smeared links can have greatly improved chiral properties. In a 

recent paper we introduced a simple and effective method to simulate four flavors of 
staggered smeared link fermions. In this work we extend the four flavor method to two 
flavors. We define the two flavor action by the square root of the four flavor fermionic 
determinant and show that by using a polynomial approximation the two flavor action can 



be evaluated with the necessary accuracy. We test this method by studying the finite 



temperature phase structure with hypercubic smeared (HYP) link staggered action on 



^ ! N t = 4 temporal lattices. 
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I. INTRODUCTION 



In a series of recent papers we introduced the hypercubic blocked (HYP) action where the 
gauge links that connect the fermions are smeared with an optimized local block transfor- 
mation, the hypercubic blocking [1]. We showed that even one level of hypercubic blocking 
improves the flavor symmetry of staggered fermions by an order of magnitude but it distorts 
the local properties of the configurations only minimally. The static potential calculated 
with HYP blocked links, for example, is indistinguishable from the thin link one at distances 
r/a > 2. We proposed a simple but effective method to simulate four-flavor dynamical stag- 
gered actions with smeared links and presented the first results on the finite temperature 
phase structure with HYP action in Ref. [2]. 

Our thermodynamic results with the HYP action are very different from the well known 
results obtained with thin link four flavor staggered fermions but consistent with the pre- 
dictions of the instanton liquid model [3]. The thin link action shows a pronounced first 
order phase transition on N t = 4 temporal lattices. Two state signals have been observed 
also on N t = 6 and 8 lattices though the discontinuity weakens as the temporal lattice size 
increases. The critical temperature is estimated to be T c ~ 150 — 170 MeV [4]. In contrast 
to that with the the HYP action we did not find any sign of a first order phase transition 
neither on N t = 4 nor on N t = 6 temporal lattices even with very light quarks. The quark 
mass dependence of the chiral condensate suggests that if a phase transition occurs at all 
in the chiral massless limit it is at a much lower temperature value than the above quoted 
T c and simulations would require larger values of N t even with the HYP action. We argued 
that the difference is due to the improved flavor symmetry of the HYP action and suggested 
that the thin link phase transition with N t = 4 — 8 is not much more than a lattice artifact. 

The simulation method for smeared actions presented and used in Ref. [2] works with any 
fermionic formulation but requires the (stochastic) evaluation of the fermionic action. That 
can be done straightforwardly for four flavors of staggered or two flavors of Wilson/clover 
fermions but requires further considerations for one or two flavors of staggered fermions. In 
this work we approximate the two flavor staggered fermion determinant with the square root 
of the four flavor fermionic determinant and evaluate the square root by using a polynomial 
approximation. Because the fermions couple to the smooth HYP fat links even a low order 
(32-64) polynomial is sufficient to evaluate the fermionic action accurately. We use this 
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simulation technique in a preliminary study of the finite temperature phase structure of the 
two flavor staggered HYP fermions on lattices with N t = 4 temporal extension. 

Before describing our method to simulate the two-flavor HYP action we would like to 
mention an issue general to any two and 2+1 flavor simulations based on staggered fermions. 
The partition function of a dynamical system with rif number of fermion flavors is 

Z nf = J DUDipDip e -W>-£rii rn ( v,m )r (1) 

= J DUdet{Q n f{V,m))e- SG ^ (2) 

where Sq{U) is the gauge action depending on the thin gauge links C/j i(U and Q(V,m) is the 
fermion matrix for one fermionic flavor depending on the smeared links V and quark mass 
m. The dependence of the smeared links on the thin links is arbitrary but deterministic. 
For rif — 4 flavors of staggered fermions the fermionic matrix is 

'.//;) = MM 

where 

M(V, m)ij = 2m5 itj + ]T ViA^A,^ ~ ^, M <W ■ ( 3 ) 

Here r] itlM are the usual staggered fermion phases and M^M is defined on even or odd sites 
only. The naive generalization for arbitrary flavors suggests the fermionic matrix 

n n t(V,m) = (M f M)^ /4 (4) 

and the partition function 

Z nf = J DU det((M^M) n f^) e- SG{u) (5) 

= J DUD^Dije- 3 ^-^ . ( 6 ) 

While eq. (5) is formally well defined, its physical meaning is less obvious for rif that is 
not the multiple of 4. Even if M^M described four degenerate flavors, its fractional power 
is a non-local quantity and the action in eq. (6) is non-local. To make the matter even 
worse, the staggered fermion action does not describe four degenerate flavors. It has only 
a remnant U(l) flavor symmetry and only one Goldstone pion. The symmetry structure of 
the fractional power of the staggered fermion matrix is not at all clear. Nevertheless actions 



3 



like the one described in eqs. (5-6) have been used extensively to study two and 2+1 flavor 
systems. One might argue that in the continuum limit there must exist a local action that 
describes a single fermion species and the fourth power of this action should agree with 
the four flavor staggered action. In that case the l/4th power of the staggered fermion 
matrix can be considered as a non-local approximation to the well defined one-flavor action. 
Perturbation theory supports this argument [5]. At finite lattice spacing non-perturbative 
lattice effects can destroy this approximate agreement and simulations with rif staggered 
flavors cannot be considered as true rif flavor simulations. How can one test if and to what 
accuracy eq. (5-6) describe rif degenerate dynamical species? There are two sides of this 
question. First, even with rif — 4 when the lattice fermion action is local, we have to ask if 
flavor symmetry is sufficiently restored so the staggered action can be considered to describe 
four degenerate flavors. Second we have to consider when it is a good approximation to 
take a fractional power of the determinant to describe rif = 1 or 2 flavors. One does not 
expect the rif ^ 4 simulations to be correct when the corresponding rif — 4 simulations are 
not, but there is no guarantee that lattice artifacts do not increase significantly when the 
fractional power of the determinant is taken. 

The properties of the vacuum give the cleanest signal to answer both questions. For 
example the topological susceptibility at finite, small quark mass is inversely proportional 
the the fermion number and the proportionality constant is predicted by chiral perturbation 
theory. Presently available data suggests that with thin link staggered fermions the lattice 
has to be fairly smooth, a < O.lfm to describe two flavors of degenerate fermions. On lattices 
with a ~ 0.1 7fm both the two and four flavor simulations fail at reproducing the correct 
chiral behavior. With the HYP action already at a ps 0.1 7fm the topological susceptibility 
is close to its chiral value, at least with four flavors [6]. 

We will not be concerned with the above outlined theoretical problem concerning two 
flavor staggered fermions in this paper. We consider the action eqs. (5-6) and describe 
an effective and simple method to simulate it. Whether this action describes the correct 
number of flavors has to be investigated for different actions at different lattice spacings 
independently. 
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II. SIMULATING TWO FLAVORS OF STAGGERED FERMIONS 



In Ref. [2] we simulated actions of the form of eq. (2) using a two step algorithm. First 
a subset of the thin links are updated using an over relaxation or heat bath updating based 
on the pure gauge action Sq{U). The fermionic determinant is not invariant under such an 
update and the proposed change is accepted only with probability 

Pacc(V',V) = mm{l,exp[-S F (V) + S F (V)]} , (7) 

where V denotes the new fat link configuration and S F (V) = —trln(Q nf )) is the fermionic 
action. If a given sequence of thin link updates is generated with the same probability as 
the reverse sequence this procedure satisfies the detailed balance condition. The fermionic 
action can be evaluated using a stochastic estimator if the fermionic matrix has the form 

SFf = Q\V)Q(V) . (8) 

In that case the acceptance probability eq. (7) can be approximated stochastically as 

PL(V, V) = min {l, exp [Q\V')Q{V) - Q\V)Q{V)] t)} , (9) 

where the vector £ is generated according to the probability distribution 

P(Oocexp{-£tgt ( y/ )g( y^} . (10 ) 

For rif — 4 flavors of staggered fermions the algorithm is straightforward. With the sub- 
stitution Q — M, where M is the matrix defined in eq. (3), the vector £ can be written 

as 

£ = M~ l R , 

where R is a Gaussian random vector and the action difference is easily evaluated as 

-S F {V) + S F (V) = ^ [M\V)M(V') - M\V)M(V)\ i . 

In order to simulate rif = 1 or 2 flavors we have to write Q nf = (M^M) n ^/ 4 in the form 
of eq. (8). That can be done using a polynomial approximation for fractional powers. In 
case of rif = 2 we need 

x 1 ' 2 = P 1/2 (x) = lim^P^lix) (11) 
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or 

,1/2 = _ p , ,J„\ = „Um 



x / = iP. 1/2 (x) = xZim n _ >00 Pli /2 (x) , (12) 

/2 V 



where -P^W^) = X^Lo^™^ i s an n th order polynomial approximation of the function 



x ±i/2 rp^g f orm f e q ^2) is more convenient to use as the polynomial P^/ 2 ( x ) can ^ e 
chosen such that it has only complex roots if n is even and can be written as 

n/2 
i=i 

where denotes the ith complex root with positive imaginary part and 

?%(x) = V^n%-rf B) ). (13) 

i=l 

With this notation the fermion matrix for rif — 2 staggered fermion flavors is 

<> 2 = (M f M) 1/2 = M t MP_i/ 2 (M t M) 

= ^m^g^MtM) M^Mq { ^ /2 (M^M) , 



-1/2 



where we have used the fact that q^ /2 (M^M) commutes with M^M. Identifying 



Q = lim n ^ 00 Mq i ^ /2 (M^M) (14) 

in eq. (8) allows us to simulate the two flavor staggered action. The vector £ of eq. (10) 
can be generated as 

£ = Umn^q^iM* M')R , (15) 

where R is a Gaussian random vector and the action difference in eq. (9) can be calculated 

as 

AS = -S F {V) + S F (V) 
= Um n ^ [P^ /2 (M'^M') M*M' - P^ /2 {M^M) M f m] £ (16) 

III. IMPLEMENTATION OF THE POLYNOMIAL APPROXIMATION 

We followed the method outlined in Ref. [7] to construct the polynomial approximation 
for the function x a in the form P^\x) = Y%=o c\ n ^x l . To find the coefficients cf^we minimize 
the integral 

/= / dx{x a - P {n \x)fx w . 
Jo 



The term x u regularizes the x ~ behavior of the integral and has to be chosen such that 
2a + uj > 0. With the limits of the integral set to and A the resulting polynomial will 
approximate x a in the (0,A) interval. The minimization condition of the integral I leads to 
n + 1 linear equations that can be solved analytically. We used Mathematica to determine 
the coefficients cf 1 ^ for a — — 1/2 at several u values for polynomials of order n = 32 — 256. 
The roots of the polynomials can also be determined using Mathematica. The roots r and 
the coefficient for several n values and uj = 1 can be found on the web site http:/ /www- 
hep. colorado.edu/~anna/Polynomials/. The coefficients on the web site correspond to A = 1 
and have to be rescaled if a different A is required. 

In staggered fermion simulations we need to approximate (M^M)~ 1//2 in the interval (0, A) 
that covers the eigenvalue spectrum of the matrix M^M. For free staggered fermions the 
maximum eigenvalue is A = 16. In simulations we used A = 18 to cover the occasional 
A > 16 eigenvalues that arise due to fluctuations. 

In evaluating £ or AS we need to multiply lattice vectors by q^\/ 2 or q^J 2 - Each such 
operation involves n/2 multiplications by M^M. A typical example is 

, n/2 

$ = q%{M^M)cp = Vc£° Hi A/ A/ - rfV (17) 

i=i 

where f is an arbitrary source vector and $ is the resulting vector of the operation. To 
reduce numerical round-off errors it is important to order the terms in the product such 
that the norms of the intermediate lattice vectors do not fluctuate too much [8]. One such 
ordering is also given on the above web site. It is also helpful to factor out the scale A and 
rewrite eq. (17) as 

n/2 A/ft M r (n) 

i'=l 

where i' denotes some chosen sequence of the orig inal index i and = (c< b) ) 1/b . This way 
the multiplication with q^y 2 can ^ e performed without excessive round-off errors even with 
n = 256. 

To simplify the notation in the following we will not write out the argument M^M or 
M'W in the polynomials but use the notation g_ 1/2 = g_i /2 (M t M), q'_ 1/2 = g_i /2 (M' t M / ), 
etc. 

By using a finite order polynomial approximation of q-1/2 we introduce systematical 
errors. One can reduce these errors both by increasing the order of the approximating 
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polynomial and by improving the operator used in evaluating AS in eq. (16). To see the 
latter let us denote the action difference evaluated according to eq. (16) using an nth order 
polynomial by AS^ and use the notation = Pi"l 2 — P-i/2- Now we can write the 
action difference as 



AS (n) = AS + ^ [A'W M*M' - A^ M f M] £ . 



(18) 



Since AS is small (otherwise the proposed change is not accepted), M'^M' — M^M will also 
be small. That implies that the correction to AS is small if A^ is small. But we can do 
even better by considering the quantity 

AS (n) = e \{P%Y (M'tM') 2 - (P%f (MtM) 21 

= AS + 3ft [A'W M' { .\l' - A'"' M*m] £ + 0(A^ 2 ] 



The combination 



AS™ = i(3AS (n) - AS (n) ) = AS + 0(A™ 2 ) 



has only 0(A^ 2 ) errors. A similar improvement is possible for the vector £ but it is too 
cumbersome and we have decided to use a higher order polynomial with the simple form of 
eq. (15) instead. 



IV. SUMMARY OF THE N F = 2 FLAVOR UPDATING 

In this section we summarize the accept-reject step of our algorithm that follows the 
sequence of over relaxation or heat bath updates of the gauge links. The over relaxation 
and heat bath steps can be performed the same way for two flavors as for four flavors and 
the details are discussed in Ref. [2]. 

Once the thin links are updated and the smeared V links are calculated the proposed 
change is accepted with a probability given by eq. (7). The acceptance probability is 
evaluated using a stochastic estimator and in Ref. [9] we showed that the acceptance rate 
can be increased significantly if the most ultraviolet part of the fermion determinant is 
removed. We define a reduced fermion matrix as 

M(V) = M r (Y)A(V) with (19) 

A(V) = exp \a A D 4 (V) + a 2 D 2 (V)] , (20) 
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where D is the kinetic part of the fermion matrix and the parameters 0:4 and a 2 are arbitrary 
but real. This way we achieve that an effective gauge action 

Seff(V) = -^ ai Retr{D 4 (V)] - ^-a 2 Retr[D 2 (V)] (21) 

is removed from the determinant. The action difference now has two parts 

AS = S eff (V) - S eff (V) + S Fr (V) - S Fr (V'), (22) 

where S Fr is the fermionic action calculated with the reduced fermionic matrix. For two 
flavors using an nth order polynomial approximation the reduced fermionic action for the 
original smeared fields is calculated as 

S Fr (V) = $ f X (23) 

with 

$ = A- 1 / 2 (V)q%(M\V)M(V))Z, 
X = M\V)M{V)i, 

or with an improved operator as 

S Fr {y) = h0X-X^X) (24) 

with 

X = P^ /2 (M\V)M(V))X . 
Here the vector £ is generated as 

£ = A(V') 1/2 q { ^f 2 {M\V')M(V')) R, (25) 

where R is a Gaussian random vector. Note that in computing the vector £ we use a 
polynomial of order m that can be different from the order n used in computing the vectors 
$, X, X . For completeness in eqs. (23-25) we again stated the argument of the function 
In deriving these formulas we used the fact that the matrices A and M}M r are 
functions of the matrix M^M and commute. The vectors £, $, X, X are defined on the even 
sites of the lattice only. The reduced fermionic action S Fr (V) for the updated smeared fields 
is calculated similarly and the proposed update is accepted or rejected with the probability 
eq. (7). The parameters a 2 and in eq. (20) can be optimized to maximize the acceptance 
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rate. We keep the same choice a-i = —0.18 and 04 = —0.006 as in [9]. It is worth emphasizing 
that this algorithm works with smeared links only. Smearing removes most of the ultraviolet 
fluctuations and the acceptance probability of the heat bath and over relaxed algorithms is 
large enough to make the algorithm efficient. 

V. THE PERFORMANCE OF THE ALGORITHM 

We have tested the above algorithm on 8 3 x 4 lattices around the critical point of the 
N t — 4 phase transition with quark masses am q = 0.01 and am q = 0.04. We found that the 
two flavor simulation is considerably more efficient than the four flavor simulation. One can 
update many more links (about twice as many) in the over relaxation and heat bath steps 
without decreasing the acceptance rate. The effectiveness of the algorithm did not change 
as we lowered the quark mass. 

To test the accuracy of the polynomial approximation we calculated the fermionic action 
eq. (23) and its improved version eq. (24) using the same £ vector and compared the results 
obtained with n = 32, 64, 128 and 256 order polynomials. Using the improved form eq. (24) 
we found very little difference between the polynomials, even n = 32 predicts the action 
difference to an accuracy of 0.01 and the acceptance probability to an accuracy of better 
than 10~ 3 on 8 3 4 lattices. The simpler form of eq. (23) gives similar results with n = 64 or 
higher order polynomials. On larger lattices higher order polynomials might be necessary. 
On an 8 3 x 24, (5 = 5.3, am q = 0.04 lattice that has a lattice spacing a m 0.23fm and physical 
volume over to 30fm 4 we found it necessary to use 64th order polynomials even with the 
improved form. 

To check the accuracy of the £ vector we performed two long runs using m = 64 and 
m = 128 order polynomials in eq. (25) at a coupling deep in the confining region on N t = 4 
lattices (f3 = 5.0, am = 0.01). We found no difference between the two runs within statistical 
accuracy. 

VI. FINITE TEMPERATURE SIMULATIONS 

We tested the algorithm presented in this article by studying the finite temperature phase 
diagram of two-flavor QCD. Our first results were obtained on 8 3 x 4 lattices for two values 
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Figure 1: On the left, the chiral condensate < jpip > divided by the bare quark mass is shown as 
a function of j3. The results are from simulations on 8 3 x 4 lattices at two quark mass values. On 
the right, E, the condensate extrapolated to zero quark mass according to eq. (27) is plotted. 



of the bare quark mass, am q = 0.04 and am q = 0.01. We consider the chiral condensate 

1 



trM~\V) , 



(26) 



where N s is the spatial lattice size. The chiral condensate divided by the bare quark mass 
is shown in the left plot of Fig. 1. The two mass values agree for j3 > 5.2 but deviate at 
smaller couplings. At small quark masses one expects that the chiral condensate depends 
linearly on the quark mass 



+ c 



(27) 



where £ is the value of the condensate extrapolated to zero quark mass. Such an extrap- 
olation is meaningful only if all extrapolated points and am q = are in the same phase. 
In that case S is zero in the chirally symmetric phase and finite, positive in the chirally 
broken one. In the right plot of figure 1 we show E from our data. Since we have done 
simulations with only two quark mass values we have no control over the quality of the 
extrapolation nor can we be certain that our data points are in the same phase as am q = 0. 
The plot nevertheless suggests that for (3 > 5.2 we are in the chirally symmetric phase while 
for (3 < 5.2 at least one of the mass values are in the chirally broken phase. A more careful 
analysis with several quark mass values must be done to resolve the phase diagram as has 
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been proposed in Ref. [2]. Nevertheless even with this exploratory study it is evident that 
the two and four flavor systems are very different. Further work to determine the phase 
diagram and the corresponding scale of the two flavor system is under way. 
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